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ABSTRACT 

We report the first three-dimensional radiation magnetohydrodynamic (RMHD) simulations of pro- 
tostellar collapse with and without Ohmic dissipation. We take into account many physical processes 
■ required to study star formation processes, including a realistic equation of state. We follow the evo- 

lution from molecular cloud cores until protostellar cores are formed with sufficiently high resolutions 
, without introducing a sink particle. The physical processes involved in the simulations and adopted 

i!^ ■ numerical methods are described in detail. We can calculate only about one year after the formation of 

the protostellar cores with our direct three-dimensional RMHD simulations because of the extremely 
QQ \ short timescale in the deep interior of the formed protostellar cores, but successfully describe the 

I early phase of star formation processes. The thermal evolution and the structure of the first and 

second (protostellar) cores are consistent with previous one-dimensional simulations using full radi- 
ation transfer, but differ considerably from preceding multi-dimensional studies with the barotropic 
approximation. The protostellar cores evolve virtually spherically symmetric in the ideal MHD models 
' because of efficient angular momentum transport by magnetic fields, but Ohmic dissipation enables 

the formation of the circumstellar disks in the vicinity of the protostellar cores as in previous MHD 
studies with the barotropic approximation. The formed disks are still small (less than 0.35 AU) be- 
cause we simulate only the earliest evolution. We also confirm that two different types of outflows are 
^ I naturally launched by magnetic fields from the first cores and protostellar cores in the resistive MHD 

\-{ , models. 

C/j ■ Keywords: ISM: clouds — ISM: jets and outfiows — magnetohydrodynamics (MHD) — radiative 

' transfer — stars: formation 
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I. INTRODUCTION 



> 

— . As stars are the most fundamental elements in the universe, their formation is one of the longstanding key topics in 
\^ ' astrophysics. Because of its coniplicate d nature and observational difficulties, computational simulations have played 
' crucial roles in this field. iLarsonl (|1969 f) first showed the scenario from molecular cloud cores to protostellar cores with 
one dimensional hydrodynamic simulations using the diffusion approximation for radiation transfer, and found that two 
' remarkable quasi-hydrostatic objects are formed in the course of protostellar collapse. A molecular cloud core initially 
\ collapses almost isothermally because thermal emission from dust grains is efficient in this phase. The gas temperature 
. starts to rise when the central region gets dense and radiation cooling becomes inefficient, then the collapse almost 
■ stops and a quasi-hydrostatic object forms. This object, a so-called first (hydrostatic) core, evolves through accretion 
from the envelope. When the central temperature reaches about 2,000 K where hydrogen molecules start to dissociate, 
the core becomes unstable and collapses again because H2 dissociation is strongly endothermic and its gas pressure fails 
to balance gravity. This second collapse ends when molecular hydrogen dissociates completely and another quasi-static 
object is formed. This is the second core or the protostellar core. It acquires its mass via accretion on the dynamical 
, timescale of the natal cloud core an d evolves into a protostar an d eventually a main-seq u ence s t ar. More sophisticated 
■ ■ simulations have been performed (jWinkler fc NewmanI Il980al lbl: iStahler et al.l llOSOal fbl. 119811 : iMasunaga fc Inutsukal 
[2001 and the scenario of protostellar collapse under spherical symmetry is by now well established. 

However, there is non- negligible internal motion (i.e., rotation and turbulence) in molecular cloud cores which makes 
the collapse completely different from the spherically-symmetric models. Because centrifugal force is typically strong 
enough to prevent formation of stars ("the angular momentum problem"), there must be sufficiently effective mech- 
anisms of angular momentum transport. Gravi t ationa l torque via non-axisymmetric structure like spiral arms takes 
this role if magnetic fields are very weak. iBatd ((loos') first performed 3D Smoothed Particle Hydrodynamics (SPH) 
simulations and showed that gravitational torque re distrib utes angular momentum sufficiently and enables formation 
of protostellar cores. On the other hand, iTomisakal (jl998l 120001 12002D showed that magnetic fields transport angu- 



lar momentu m far more efficiently via torsional Alfven waves (so-called magneti c braking. iMouschovias &: Paleologoul 
()1979l 119801) ) and bipolar outflows using two-dimensional MHD simulations. iMachida et al.l ( 2005bl laf) performed 
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three-dimensional nested-grid MHD simulations of simila r situations and invest igated the condition for fragmenta- 
tion in the early phase before the second collapse begins. iMachida et alj (|2006aD investigated the evolution after the 
second collapse and showed that two different outflows are launched from the first core an d the protostellar core. 
More detailed 3D simulations on this problem have been performed (Ba neriee fc Pudr itz 2006; ^ Hennebelle fc Fromangl 
[2OOI IHennebelle & Tcvssicr''2008'; 'Machida et al."2008b') and effects of many physical ingredients s uch as m isaligned 
magnetic fiel ds (Matsu inoto &: Tomisaka 200 4: . Machida et al.. ,2006b: ,Ciardi &: Hcnn cbcUe 201ft iJoos et al. 2012). 
metallicities ()Machidall2008f ) and turbulence~ " jSeifried et al. r i2012l: ISantos -Lima et al.ll20i2f ). have been investigated. 
Even weak magnetic fields can transport angular momentum very efRciently and modify the stru cture of the accre- 
tion flow significaiitly, and observations suggest that molecular clouds are considerably magnetiz ed (jHeiles fc Crutcheil 
l2005t iGirart et al.l[2Ci06l : iFalgarone erani200l iTroland fc Crutchedl200l iCrutcher etHllBol . Therefore, magnetic 
fields are likely to play dominant roles in typical present-day environments. 

Non-ideal MHD effects are also important in protostellar collapse. Typical magnetic flux in a molecular cloud core 
is far larger compared to that i n a formed star, which me ans that there must be significant redistribution of magnetic 
flux during collapse. Moreover, iMellon fc Lil (|2008l . 120091) claimed that rotationally-supported circumstellar disks are 
difficult to form in magnetized molecular cloud cores because the angular momen tum transport due to magnetic 
fields is too efficient ("the magnetic braking catastrophe", see also iLi et al.l ()2011[ )). Non- ideal MHD effects such 
as Ohmic dissipation and ambipolar diffusion are supposed to be responsible in this process, because the ionization 
degree in star forming clouds is very low. Ohmic dissipation redistributes sufficient magnetic flux outward from 
the centrally-conden s ed high density region and enables format i on of circumstellar disks in th e vicinity of protostars 
papp fc Basul 120101: Unutsuka et aT]|2010t [Machida et al.ll2011t IMachida fc Matsumotoll2011[ ). Gravi tational torque 
becom es important in transporting angular momentum where magnetic fields are significantly weakened ([Machida et al.l 
[20Tnl) . 

In those previous multi-dimensional simulations, they adopted the barotropic approximation for gas thermody- 
namics in which the the rmal evolution is approximat ed with simple polytropic relations mimicking the results of IT) 
RHD simulations (e.g., 'M asunaga fc Inutsukal [20001) . instead of solving radiation transfer. Recently, iBatd (|2010l 
see also [Wh itchous efc Batd 120061 ) reported 3D SPH RHD simulations of formation of protostellar cores. 
Schonke fc Tscharnuten (|20l"ll) also performed similar calculations of protostellar collapse using 2D axisymmetric RHD 
simulations. They commonly found interesting phenomena happen when the protostellar cores are formed; bipolar 
outflows are launched from the disk-like flrst cores via radiation heating (not radiation force) from the protostellar 
cores. Considering the energy released in the second collapse and subsequent accreti on onto the protostellar core, the 
irradiation is sufficient to heat up the gas in the first core and launch the outflows ([Batel[20TT[ ). The barotropic ap- 
proximation cannot reproduce such outflows driven by radiation heating. The structure of the flrst core which depends 
on the angular momentu m distribution seems to be im portant in this phenomenon, and in the extreme case, the flrst 
core is almost blown up (Schon ke fc Tscharnuterl[2011[ ). Such violent phenomena may affect the story of protostellar 
collapse, circumstellar disk formation, and evolution of protostars. 

Thu s bot h magnetic fiel d s and r adiation transfer are critically important in star formation processes. iCommergon et all 
([2OIOD and iTomida et al.l ([2010bD independently performed 3D RMHD simulations of protostellar collapse, but both 
stopped their calculations before the onset of the second collapse and they adopted the ideal MHD approximation. In 
this work, we report the flrst 3D RMHD simulations of protostellar collapse from molecular cloud cores to protostellar 
cores with and without Ohmic dissipation. To follow the evolution after the second collapse, we adopt a realistic 
equation of state (EOS) which takes chemical reactions into account. The goal of this study is to describe the realistic 
evolution in the early phase of protostellar collapse (i.e., until protostellar cores are formed) involving both radiation 
transfer and magnetic flelds. 

The plan of this paper is as follows: We describe the basic equations in Section 2. The adopted numerical methods 
are given in Section 3 and the simulation setups in Section 4. We show the results of our 3D RMHD simulations 
in Section 5. Section 6 is devoted to conclusions and discussions. In the Appendices, we explain the microphysical 
processes considered in our simulations. 



2. BASIC EQUATIONS 

In this work, we solve system equations including MHD, self-gravity, radiation transfer, a realistic EOS, and Ohmic 
dissipation. We use the comoving-frame RHD equ ations in which all the radiation related variables are deflned in the 
comoving frame of the local fluid element (Castori[2003). For simplicity and to reduce computational load, we adopt 
the gray approximation (using frequency-a veraged equations assuming the local thermodyn amic equilibrium) and the 
flux limited diffusion approximation (FLD, iLevermore fc Pomranindl 19811 : iLevermoril 1 984[ ) for radiation transfer. We 
take Ohmic dissipation into account as t he first step because it is supposed to be the most significant in the context of 
magnetic flux loss from high density gas (IMachida et al.ll2007D. alt hough other non-ideal MHD effects such as ambipolar 
diffusion and Hall effects are also important ([Nakano et al.l[2002l ). 

The governing equations are as follows: 



dp 



dpv 



pv ( 



B(g)B 



V • (pv) = 0, 



c 



(1) 

(2) 
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— - V X (v X B - r)V X B) = 0, 
dt ^ 11^ 

V • B = 0, 



de 

dt 



dEr 



e+p+ -|Bp ) V - B(v • B) - 7/B x (V x B) 



V • [wEr] + V • F^ + P,. : Vv = cap{arT^ - E,.), 



(3) 
(4) 

(5) 
(6) 

(7) 



}Er 



R 



WeA' 



where p denotes the gas density, v the fluid velocity, B the magnetic flux density, p the gas pressure, Tg the gc 
temperature, e = Cg + ■^pv^ + ^|Bp the total gas energy density (e^ is the internal energy density of the gas 
Er the radiation energy density, $ the gravitational potential, F^ the radiation energy flux, and 



the radiation 



2.99792 X IQi^cms 



is the speed of light and G 
-"erg cm~ 



6.673 X 10- 



7.5657 X 10 ^^erg cm ^ K is the radiation (density) constant where a 



pressure tensor, respectively, c 
the gravitational constant, ar — 4(t/c 
5.6704 X 10~^ergcm~^ K""* is the Stefan-Boltzmann constant, ry is the resistivity and an (up) is the Rosseland (Planck) 
mean opacity. I denotes the unit matrix. ":" means the double dot product of two tensors, A : B = AijBji. 
These equations represent conservation of mass, the equation of motion, the induction equation including Ohmic 
dissipation, the solenoidal constraint, the Poisson's equation of gravity, and the gray FLD radiation transfer equations 
from top to bottom. Basically we use the Gaussian cgs units but we rescale the magnetic flux density to eliminate 
the constant coefficients, i.e., B = Bo/-\/47r where Bq is given in Gauss. Additionally, we need an EOS which gives 
the relations between the thermodynamic variables p, p, T and eg to close the system. We adopt the tabulated EOS in 
which we consider internal degrees of freedom and chemical reactions of seven species (H2, H, H"*", He, He"^, He^+ and 
e~). Here we assume the solar abundance, X = 0.7 and Y = 0.28. We also adopt the tabulated tables o f the Rosseland 
and P l anck mean opacities. In o rder to cover the huge dynamic range, we combine thre e opacity tables: iSemenov et ahl 
(|2003l ). lSeaton et all (|1994[ ) and lFerguson et al.l (|2005D . The resistivity is adopted fromlUmeb avashi fc Nakanol (|2009f ). 
while we additionally consider thermal ionization of potassium (K). The detailed descriptions about these microphysical 
processes ar e given in Appe i idices. 

Note that lMachida et all ()2006al [20071 l2008b[) simplified the resistive term in eq. [3]to TyV^B, but it is inadequate 
because the resistivity r] has large spatial gradient according to the local gas density and temperature gradients, and 
what is worse, it violates the divergence free condition (eq. 14]) . In this work we correctly discretize the original equation, 
and also includes the effect of the resistivity in the energy equation (eq. [5]). 

3. METHOD 
3.1. Operator-Splitting 

In order to solve the complex system described in the previous section, we divide the system into five parts and solve 
them separately on the nested-grid hierarchy. We update our system in the following strategy: 
Step 1. Ideal MHD part: 



dp 
dt 



dpv 



pv I 



2' ' 
dB 
'dt ~ 



- V • (pv) = 

:-B(g)B 

' V X (v X B): 

V B: 



de „ 



v-B(v B) 



0, 


(8) 


0, 


(9) 


0, 


(10) 


0, 


(11) 


0, 


(12) 


0. 


(13) 



Step 2. Self-Gravity part: 



V^<^ = 4:nGp, 



(14) 
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dt 

Step 3. Resistivity part: 

dB 

'dt 

de 
dt 

Step 4. Radiation part: 



^ = -pvV$, (16) 



Vx(77VxB)=0, (17) 
V • [?7B X (V X B)]=0. (18) 



dt 

Step 5. Radiation Force part: 



-^ = -capiarT^~Er), (19) 

dE 

V -Fr + Fr ■■y^ = cap{arT^ - Er), (20) 

%^ = ^F., (21) 

"«F..v. (22) 



9i 



c 



Introducing such an operator-splitting technique causes loss of time-accuracy, and the overall accuracy of our time- 
integration scheme is first order. 

3.2. Magnetohydrodynamics and Self-Gravity 

We solve the MHD system using the HLLD approximate Riemann solver (jMivoshi fc Kusaiioll2005| ). It is suitable 
for simulations involving the realistic EOS because this scheme does not require the detailed knowledge about the 
EOS, unlike Roe's approximate Riemann solver. It is simple, robust, highly efficient, and almost as accurate as Roe's 
solver. We define all the variables at the cell center. To numerically satisfy the solenoidal constraint (|4] ), wc adopt 
the mixed correction based on the generalized Lagrangian multiplier approach proposed by iDedner et al.l |2002.) . To 
achieve second-order accuracy in space and time, we adopt standard MUSCL (Monotone Upstream-ce ntered Scheni e 
for Conservation Laws) approach and the directionally-unsplit two-step predictor-corrector scheme (e.g.. lHirschlll99Cil ). 
We store (p, v, B, Cg, tp^ Er) as primitive variables and interpolate them for spatial reconstruction with the minmod 
slope limiter for robustness {ip is an additional variable related to the solenoidal constraint). We use the gas internal 
energy eg as the primitive variable instead of the gas pressure p which is the textbook notation in order to minimize 
the numerical error when we use the tabulated EOS. That is, the discretization error of the EOS table will be directly 
reflected in the results if we perform mutual conversions between p and eg (i.e., p ^ eg ^ p). We can avoid this error 
if we use one-way conversion (eg — p) only; the discretization error only affects the flux and its effect is small. 

Some high resolution (magneto-) hydrodynamic solvers show strange oscillations at strong shocks aligned to grid 
structure. This problem called the Carbuncle ph enomenon is obviously un physical. To suppress this unphysical 
oscillation, we use the HLLD— flux devel oped bv iMivoshi k, Kusanol (|2Q07f ) in the cells which potentially contain 
shocks. We adopt the method proposed in iHanawa et al.l (|2008| ) to identify such cells. HLLD— is a modified version 
of HLLD, in which the tangential velocities {vy, Vz) in the Riemann fan are replaced by the HLL averages while other 
variables are kept the same as HLLD. Sufficient but not too large additional shear viscosity is introduced by this 
procedure and it stabilizes the Carbuncle phenomenon. 

We update the advective flux of radiation energy separately using Roe's upwind method. When the radiation 
energy dominates the gas energy, the radiation pressure affects the sonic speed of the gas and the characteristics are 
modified, but we neglect this effect. In our star formation simulations (at least in low mass cases), we expect that 
we encounter such a radiation dominant region only in the deep interior of the protostellar core in the late phase of 
protostellar collapse, and we do not expect this effect will be significant. 

To solve the Pois son's equation on the nested-grid hierarchy, wc adopt the multigrid solver developed by 
iMatsumoto k, Hana wa (2003a). The solver gives second-order accurate solution in space. We integrate and 
using obtained gravitational potential. 



3.3. Ohmic Dissipation 

Because of the low ionization degree, non-ideal MHD effects such as Ohmic dissipation, ambipolar diffusion and 
the Hall effect work during protostellar collapse. Ohmic dissipation becomes dominant in the high density re- 
gion (typically p ^ 10~^^ g c m~'^) and ambipolar diffusion dominates in the low density region (jNakano et al.l [20021 : 
iKunz fc MouschoviasI 120091 I2010D , but in this work we focus on Ohmic dissipation as the first step because it is 
supposed to be most significant in the context of the magnetic flux problem and the angular momentum problem 
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(|Machida et al.|[2011l : iDapp et al.ll2012[ ). Taking account of ambipolar diffusion will be a future work. The timescale 

of Ohmic dissipation becomes shorter than the dynamical t imescale in the high density region (p ;J 10^^° g cm^"^, 
Figure 1^5)) and significant diffusion of magnetic flux occurs (|Nakano et al.l [20021 ) . Note that these non- ideal MHD 
processes do not reduce total magnetic flux but redistribute centrally-concentrated magnetic flux outward from the 
high density region (see below, Figure [TS]) , although this process i s ofte n mentioned as " magnetic flux loss" . 

We discretize ([T7]) and ([TOl) in the same way as iMatsumotol ()2011f ). Additiona lly, we introduce th e correction 
terms proportional to V • B to improve the nature of the system as described in [Graves et alj ()2008l ). Since the 
timescale of the diffusion of magnetic fields can be far shorter than the dynamical timescale (Figure [28| , here we 
adopt the Super-Time-Stepping (STS) method proposed by Alcxiadcs et al. (1996) (for astrophysical applications, see 
fO^Sullivan fc DownesI (|2006l ) : IChoi et alJ (j2009l ): iCommercon et alJ (|2011bf )) when the timestep for the resistivity part 
is shorter than the hydrodynamic timestep. STS is constructed based on the simple explicit scheme, but we can relax 
the Courant-Friedrich-Levy (CFL) stability condition by using the following sequence of sub-timesteps: 



At 



exp 



{—l + l^) COS 



2j 



N 



1 TT 

~2 



^ +1 



The longer time integration of ATsts is achieved after N sub-timesteps. 



(23) 



AT. 



STS 



N 



N 



(i + V^)2^-(i-V^) 



2N 



(1 + v^)2^ + {i-V^) 



2N 



(24) 



where At^xp is the explicit timestep which satisfies the CFL condition and is a small positive parameter which 
controls the stability and efficiency of the scheme. Smaller v gives better acceleration but STS may give unphysical 
results when it is too small. The optimal choice of v depends on the problem, but typically v ^ 0.01 seems to be 
good for parabolic problems. From (j24p . the maximum acceleration for given v compared to the explicit scheme 
{ATexp — NAtexp) is 0.5/v^ and the most efficient calculation is realized around N ~ 0.5 /y/i'. In this work, we 
adopt z/ — 0.01 and = 6. We apply this scheme repeatedly until required time update is achieved. Even in the most 
time-consuming situation (i.e., the magnetic Reynolds number is very low, p ^ 10~^gcm~'^) in our production runs, 
STS keeps the computational cost lower than or similar to the total of all the other parts (including MHD, self-gravity 
and radiation transfer) at most. 



3.4. Radiation Transfer 

The timescale related to radiation can be far shorter than that related to hydrodynamics in star formation simu- 
lations. The lightspeed c = 2.99792 x 10^"cms~^ is far larger than typical fiuid velocity in star formation, which is 
on the order of (1 - 100) x lO^cms"^. Moreover, the radiation source terms are strongly non-linear. Therefore the 
differential equations of the radiation subsystem can be very stiff. It is unreasonable to calculate such a system using 
explicit schemes or even STS. Therefore we adopt an implicit time-integration scheme which is stable regardless of the 
timescale of involved physical processes. 



3.4.1. Discretization 



We adopt the first-order backward Euler method which is simple and stable. Here we discretize p9|) and (|20p in 
one-dimension (extension to multi-dimension is straightforward) as follows: 



At 



At 



— C(Jp 

1 
Ax 



Irpn+l 



cX 



Ax 



cX 



Tpn+1 



Ax 



' r.i 



(Vv), 



(25) 



(26) 



Superscripts and subscripts denote the indices of the discretized time and space, respectively. To avoid numerical 
difficulties due to strong non-linearity in t he flux limite r and the opacities and to achieve stable integration, we adopt 
the time-lagged opacities and flux limiter ()Castoiil2004[ ). This may cause loss of the time accuracy of the scheme, but 
we confirmed that it does not matter in our production runs because the timescale of the evolution of radiation fields 
is similar to the hydrodynamic timescale and therefore well resolved. We also use the time-lagged Eddington Tensor 
D", which yields the radiation energy tensor P"+^ = D"i?"+^. 

There are some ways to evaluate the opacities and the fiux limiter at the cell interfaces, i + \- Here we follow 
IHowell k. Greenoughl (|2003f ) and adopt the surface formula which gives good flux even at a sharp surface of optically 
thick material like the surface of a first core. 



(7R,i + CFR,i+l 3Ax 



(27) 
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We evaluate the flux limiter at the cell interface X{R)^_f_i using the following equations: 

l(v£;.),+i| 



R 



E. 



Er,i 



+ 1 



[VEr] 



E., 



r,i+l 



Er 



Ax 



In three dimensions, the radiation energy gradient should be replaced with 



(28) 

(29) 
(30) 



/ 



E., 



{VEr 



k+^,j,k = 



r.i+l 



Er 



Ax 

Er,i+l,j + l,k + Er/L.j + l.k ^ -E'r,i+l,j-l,A: ^ Er,ij-i,k 

Ia'z 



(31) 



3.4.2. Newton-Raphson Iterations 

To solve the non-linear system (l25t and p6|) . we perform the Newton-Raphson iterations ([Press "etaI1[2007h . In 
this method, we search for the zero-points of the residual functions /i(X). We can find the root iteratively using the 
following matrix equation based on the Taylor expansion: 



N 



dx 



(32) 



In our system, X is the vector of the gas and radiation energy densities in all the cells: 

X — (^5,1 7 ■ • • ; ^gA^ 7 ■ • • 1 Er^l , . • . , Er^i , Er^i-^l ; ■ • •) 

The residual functions are given from ([25)1 and (|26p as follows: 



f! 



eg,i + Atcap 



n+1 



n+l 



E'J 



At 



1 
'Ax 



-cap 



cA 

-e: 



Ax 

+ {n- : {Vv),} E, 



n+l 



n+l 



cAV 

^r) 



Aa 



(33) 
(34) 

(35) 



We can rewrite (|32|) explicitly 



df! 



-Se 



n+l 



■Se 



n+l 



Q n+l g,: 



dE\ 



SE. 



•n+l 



n+l r,i 



a n+l g,'i. 
"^g.i 



+ 



df! 



-5e: 



n+l 



Q^n+l r,^ 



f9 
It ' 



(36) 
(37) 



By substituting (p6|) into ()37p , we can eliminate the equation related to gas energy (iHaves et al.ll2006l ) (this procedure 
corresponds to performing partial LU decomposition analytically). Then we obtain the matrix equation: 



dfl 



dfl df! 



df! 



SE, 



n+l 



dfl 



-SE. 



n+l 



Qpn+l r,i+l 
^^r,i+^ 



dfl xj7n+l 

dKr-1 



dfl 



df^ 

^^f!-fl (38) 



de-^ 



The derivatives are given as follows: 



df 



■.dTa 



^ = l+4At ca-ariUp,, e^+i)]3^(p„ e^+i). 



9,« 
9 



de 

df! 

dE-r 

^ = -4At ca?.a.[T,(p„ e^f )]3g(p„ e^+i). 



(39) 
(40) 
(41) 
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= 1 + At 



n+l 
r,i 

dfl cAt ( A 



ca?. + Wl : (Vv), 




A 



(42) 
(43) 



Tg{p, Cg) and ^^(p, e^) are given from the tabulated EOS. Note that this Jacobi matrix is symmetric. We can obtain 
the solution by updating Er^i and Cg^i using SEr^i and Scg^i calculated from ([55)1 and ([55]) until and (5X become 
sufficiently small. As the initial guess for E"f^ and e^^^, we adopt the solution at the timestep n. In our simulations, 

we use convergence thresholds like max < 5 x 10"'' and max < 5 x 10~^. If these thresholds 

are not satisfied after many iterations, we take substeps with shorter At and try again until we obtain the converged 
solution successfully. 

3.4.3. Linear System Solver 

In three dimensional Cartesian coordinate, the Jacobi matrix in ((38|) is a very large sparse seven-diagonal matrix. 
Moreover, because of the strong nonlinearity of the system, the matrix is not necessarily diagonally dominant. Therefore 
we need an efficient and robust sparse matrix solver and extensive computational resources to solve this large (typically 
64^ = 262144 cells per grid level) linear system. We found that the combination of the BiCGStab (stabilized bi- 
conjugate gradient) solver and the incomplete LU (ILU) decomposition preconditioner without fill-in works well for 
our problem^. 

3.4.4. Radiation Force 

We simply integrate the radiation force terms in Step 5 using the obtained solution in Step 4. These terms are 
relatively small, at least in the early phase of low-mass star formation processes. 

3.5. Nested-Grid 

In order to achieve the hug e dynamic r ange required in protostellar collapse simulations, we adopt the three dimen- 
sional nested-grid techniqu e (jYorke et a l. 1993; Yorkc & Kaisig 1995; Zieglcr & Yorkc 1997; Matsumoto & Hanaw^ 
l2003btlMachida et al.l[200l . This is a simplified version of the adaptive mesh refinement (AMR) technique. Each grid 
level consists of x Ny x Nz cubic cells. The number of the cells in one direction iV» must be a power of 2. The finer 
grid is placed around the center of the coarser grid self-similarly. The size of a finer cell is half of that of a coarser cell. 
We number the levels from coarsest to finest: I = l(coarsest), 2, i(finest). 

We calculate the coarsest grid I = 1 first and then proceed to finer grid levels because we require the boundary 
conditions at the next timestep in the implicit update. We update all the grid levels in each step described in l3.1l The 
timestep is determined by the Courant-Friedrich-Levy (CFL) criterion derived for the MHD part. All the grids share 
this co mmon tii nestep and are advanced synchronously. We adopt so-called Jeans condition proposed by lTruelove et al.l 
(|1997| ) (see also lCommergon et al.l (|2008[ )^ for refinement; we generate a finer grid to resolve the minimum Jeans length 
with 16 cells. 

For the MHD and resistivity parts, we apply the standard procedures to the boundaries between the levels of different 
resolution. That is, when we update a level /, we construct the boundaries from the coarse level I — 1 using time and 
spatial interpolation. For the spatial interpolation, we adopt linear interpolation with a slope limiter to assure the 
monotonicity. After updating all the levels, we transfer the results in the overlapped regions from the finer grid to the 
coarser grid using conserved variables and recalculate the flux in the coarser level I — 1 using the obtained flux in the 
flner level I, conserving the total flux across the level boundaries. 

3.5.1. Implicit Update on Nested-Grid Hierarchy 

The implicit time integration for the radiation transport on the nested grids is more tricky. Because all the cells 
among the whole grid levels interact with each other in a single timestep in the implicit scheme, in principle we should 
integrate all the levels at the same time treating them like one huge non-uniform grid. However, it makes the Jacobian 
matrix irregular and complicated, and the computational costs can be very expensive because a huge number of cells 
are involved in the integration. Therefore we solve each grid separately as we do in the MHD part. This enables us to 
adopt computational algorithms highly optimized for regular sparse matrices. 

When we update a grid level /, we treat it as a uniform grid. The boundary conditions at the next timestep are 
constructed using linear spatial interpolation from the coarser grid which is already updated, and treated as fixed 
(Dirichlet) boundaries. In p9)) and (|20l) . we need to estimate the internal energy and temperature, rather than the 
total energy. We use the total energy for grid interaction in the MHD part, but it causes artificial thermalization of 
the internal dispersion of fluid motion and magnetic fields resulting in unphysically high entropy and temperature in 
the region overlapped with the finer grids. To avoid this, we take the temperature average of the nearest eight cells 

''' ILU type preconditioners are very robust and efficient, significantly reducing the number of iterations required in iterative solvers. 
However, it does not easily fit parallelization because of the dependencies between the operations. Therefore, although it is a good 
algorithm for supercomputers with high single node performance like NEC SX-9 which we use, its scalability is problematic for massive 
parallel architectures. 
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in the overlapped finest levels and calculate the internal energy using the gas density and the averaged temperature. 
For consistency between the gas and radiation, we also take the radiation temperature averaged over the nearest eight 
cells in the finest level. These procedures violate conservation of the total energy, but in the diffusion approximation 
it is more important to properly estimate the temperature at the cell center and its gradient, rather than the total 
energy. And its effects are kept small because we overwrite the overlapped regions with the results obtained in the 
finer levels at the end of every timestep, which satisfy the local conservation laws. 

We only consider the interaction between the coarse and fine levels at the level boundaries in our implicit scheme on 
the nested grids, but it may cause the loss of consistency and accuracy. In our production runs, we confirmed that this 
treatment only causes minor discrepancies between the grid levels in the extremely optically-thin regions. Because we 
are mainly interested in the evolution of the condensed objects, we can tolerate these errors. 



3.6. Boundary Conditions 

As the outer-most boundary conditions for the magnetohydrodynamics and radiation parts, we set all the cells 
outside the initial Bonnor-Ebcrt sphere to maintain their initial values, mimicking an isolated molecular cloud core 
confined in a static environment. For the Poisson's equation of self- g ravity, we compute the gravit ational potential 
of the isolated system at the boundaries by the multi-pole expansion (jMatsumoto fc Hanawal[2003af ). Our boundary 
conditions allow the gas to inflow into the computational domain through the boundaries, but the inflowing mass 
during the simulation is sufficiently smaller than that of the total mass of the initial cloud, about 8% in all the models. 

4. SIMULATION SETUPS 

We use unstable Bonnor-Ebert ()BonnoHll956l : lEbertlll955l hereafter BE) spheres of T = lOKas the initial conditions 
of our simulations, mimicking isolated molecular cloud cores. We construct a critical BE sphere (with a dimension less 
radius of ^ = 6.45) and make it unstable by increasing the gas density by a factor of ^o- Then we introduce uniform 
(rigid-body) rotation, uniform magnetic fields and m = 2 density perturbation, where m is the number of longitudinal 
modes. The initial density profile is given as follows: 



P{r) = PcMr){l + Ao) 



I + A2— jcos(20) 
-ft" 



(44) 



where pc is the gas density at the center of the cloud core, poir) the normalized density profile of the critical BE 
sphere, R the radius of the critical BE sphere, A2 the amplitude of to = 2 perturbation, respectively. In order to 
minimize the effect of initi al resolution, we adopt this "regularized" m — 2 perturbation which is smooth at the center 
of the cloud in contrast to lBoss fc Bodenheimeii ()1979[ ). 

In this work, we adopt pc = 10~^^gcm~^ and Aq — 0.2, which yield an unstable BE-like sphere whose mass and 
radius are M ^ IMq and R ~ 4.25 x 10"^ pc ~ 8800 AU. The initial free fall time at the center of the cloud is 
iff - 6.08 X lO'^yrs. 

We calculated a spherically symmetric model and four magnetized rotating models with and without Ohmic dis- 
sipation; the parameters of the models are summarized in Table [TJ The first letter of the model name denotes the 
treatment of magnetic fields: / denotes an ideal MHD model and R a resistive MHD model. The second letter rep- 
resents the initial rotation speed: F is fast and S is slow. Note that we choose these parameters so that the first 
core disks do not fragment before the second collapse because of the limitation of our nested-grids. Model SP is the 
spherical model without rotation and magnetic fields. For magnetized models, we impose the uniform magnetic fields 
of 20 pG parallel to the rotation axis. The corresponding mass-to-fiux ratio normalized by the critical value of stability 

is po = (M/i') t ~ 3.8 where $ = ttR-^Bq and (M/$)crit = ihY^^- Here we adopt the critical mass-to-fiux 
ratio of iMouschovias fc Spitzeii ()1976D but we should regard this value as just a guide because our initial conditions 
are not uniform. There is another similar threshold for stability between gravity and magnetization derived for disks 

(|Nakano &: Nakamural[l97l : we have Aq = = 7.2 at the center of the cloud, where {E/B)„it = {^■k^G)-'^^'^ . 

These mass-to-flux ratios indicate that our magnetized models are in the magnetically super-critical regime but con- 
siderably magnetized. 

We define the origin of time as the epoch of formation of the protostellar core, i.e., when the central gas density 
exceeds pc = lO^'^gcm^'^ for the first time, for descriptive purpose. We stop our simulations when the central 
temperature reaches Tc ~ 10^ K. Our simulations show the formation and earliest evolution of protostars. The typical 
spatial resolution around the surface of the first cores (r ~ 3 — 5 AU) is Aa; ~ 0.14AU (L = 12) and the finest resolution 
at the end of the simulations is Aa; ^ 6.6 x 10^^ AU ^ O.O14i?0 (L — 23). The nested-grid technique enables us to 
achieve such a huge spatial dynamic range of more than eight orders of magnitude. 



5. RESULTS 



5.1. 



Spherical Model 

We first show the results of the spherical model SP to understand the whole evolution from a molecular cloud core 
to a protostellar core, and to demonstrate validity of our code. 



5.1.1. Thermal Evolution and EOS 
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Table 1 

Summary of the initial model parameters. From left to right: the normalized angular velocity, the angular velocity, the magnetic field 
strength, the averaged mass-to-flux ratio, the local mass-to-fiux ratio at the cloud center, the amplitude of m = 2 perturbation and 
whether the resistivity is introduced or not. Other parameters are common: M = 1 Mq, R ~ 8800 AU, pc = 1.2 X 10~^* g cm~^, and 

To = 10 K. See the text for details. 




Gas Density [g cm"^] 

Figure 1. The evolution track of the central gas element (white) in SP overplotted on the distribution of the adiabatic index F in the 
p — T plane. Four low-F (blue) bands correspond to the endothermic reactions of the dissociation of molecular hydrogen, the ionization of 
hydrogen, the first and second ionization of helium, from bottom to top. 

Before discussing the global structure and evolution, we explain the thermal evolution of the central gas element. 
To see the effects of the microphysics involved in the EOS on the thermal evolution, we show the evolution track of 
the central gas element in the p — T plane and the adiabatic index T in Fi gure [U The evolution is cons istent with the 
results of the more elaborate ID spherically-symmetric RHD simulations (|Masunaga fc Inutsukal 120001 ) except for the 
details of the EOS. From this figure, we can understand the relation between the thermal evolution of the gas and the 
microphysical processes. 

While the gas density is low (pc ^ 10~^^ gcm""^), the gas collapses isothermally because radiation cooling is very 
efficient, and the EOS does not matter in this regime. When the gas density gets sufficiently high and radiation 
cooling becomes inefficient, the gas evolves quasi-adiabatically and the temperature rises. Beyond this point, the EOS 
plays almost dominant roles in the thermal evolution of the gas (at the center of the cloud; note that the gas in outer 

region is afi^ected by radiation transfer). While the gas is still cold (T^ ^ 100 K), the adiabatic index F is about 5/3 
because the rotational transitions cannot be excited and molecular hydrogen behaves like monoatomic molecules. F 
decreases to ~ 7/5 when the gas becomes warm enough to excite rotation. When the temperature exceeds Tc ~ 2000 K, 
molecular hydrogen starts to dissociate and the second collapse begins. The endothermic reaction of hydrogen molecule 
dissociation significantly decreases the effective adiabatic index to F ~ 1.1, below the critical value of stability of a self- 
gravitational sphere, Fdit = 4/3. After the completion of the dissociation, the gas evolves quasi-adiabatically again, 
while the ionizations of hydrogen and helium slightly affect the evolution making the adiabatic index softer. From 
Figure [1] we can see that the ionizations of hydrogen and helium do not cause the "third collapse" ; those endothermic 
reactions proceed gradually because the reverse reactions (recombinations) occur rapidly in the high density regions 
and the adiabatic index remains larger than the critical value. 
In this work we assume the ortho:para ratio of molecular hydrogen to be 3:1, but if we adopt the equilibrium ratio, 
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Figure 2. The evolution of ttie gas density (top) and temperature (bottom) profiles in Model SP. The numbers show the order of evolution. 



additional ener gy is consumed to convert parahydroRen to ortho hydrogen {AE/k ^ 170 K), resulting in the softer 
adiabatic index ((Boley et al.ll2007l : iStamatello s fc Whitworthll2009() . The gas does not forget the thermal history in this 
phase because its evolution is almost adiabatic, therefore the ortho-para ratio has non-negligible impact on the whole 
star formation process. We have to keep it in mind that the differences in the EOS quantitatively affect the evolution 
and properties of the first and protostellar cores such as their radii, masses and stability against fragmentation. 

5.1.2. First Core 

The first core is formed at t ~ — 650yr when the central density exceeds pc ^ lO^^^gcm"'^. Its radius is initially 
about ~ 5.4 AU, and it contracts gradually to ~ 3AU as it evolves (Figure [2]). The evolution and structure of the 
first core are in good agreement with the results o f the ID spherically s ymmetric RHD simulations using full non-gray 
radiation transfer (Masunaga ct al. 1998; Masun aga &: Inutsukal [2000l ) and the 3D SPH RHD simulations using the 
gray FLD approximation f|Whitehouse fc Batell20()6h . The outer region attains a higher entropy than the central gas 
element because of heating at the shock and radiation from the central hot region (Figure [3]). The shock at the surface 
of the first core is isothermal, i. e., it is a super-cri t ical sho ck at which almost all the kinetic energy is radiated away 
to the upstream, as discussed in iCommercon et ahl (|2011aD . The lifetime of the spher ical first core is a bit longer than 
that in iBatd ()201lD . probably because of the different initial conditions; iBatd (|2011t ) adopted a uniform sphere with 
higher central gas density, which is more unstable and gives higher accretion rate than our BE sphere. 

When the gas temperature exceeds the evaporation temperature of all the dust components (T ~ 1500 K), the 
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Figure 3. The distributions of the thermodynamic properties in the p - 
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T plane. The evolution traclc of the central gas element is also 



opacities around the central region drop significantly. The temperature distribution within this dust free region 
becomes almost flat. As Figure [2] indicates, the dust evaporation front is located at i? ^ 1.2 A U at the end of the first 
core p hase (t ~ — 2.7yr). The dust evaporation front is also visible as a small jump in Figure|31 ISchonke fc Tscharnuteri 
(|2011| ) pointed out the importance of the dust evaporation on the dynamics of first cores, but it is not prominent in 
the spherically symmetric model. 

5.1.3. Protostellar Core 

The second collapse begins when the central temperature exceeds Tc ^ 2000 K. Soon after the onset of the second 
collapse, the protostellar core is formed in the short dynamical timescale of several years (the free fall time corresponding 
to the central density when the second collapse starts {pc ~ 10~*gcm~^) is only about 0.67 years). Within 0.7 years 
after the formation of the protostellar core, it acquires Mpc ~ 2 x 10~^ Mq and the averaged accretion rate is very 
high, 2.7 X lO~^M0yr~^. The protostellar core expands due to the addition of newly accreted gas and the accretion 
shock at the surface of the protostellar core quickly propagates outward. The jump condition at this shock is almost 
adiabatic, which means that the flow is radiatively inefficient ("hot accretion") in this early phase. The outer region of 
the protostellar core is heated up by the shock and attain a high entropy, leading the protostellar core to be convectively 
stable (Stabler et al. 1980a. b). These thermal properties and expansion of protostellar cores cannot be reproduced by 
the barotropic approximation in which the shock heating is not taken into account. At the end of the simulation, the 
radius of the protostellar core is about Rpc ^ 0.047 AU ^ lOi?0. From the virial theorem, the energy released in this 



phase can be estimated to be 



PC 



hydrogen, 



XM 



2B.PC 



6.8 X 10 erg, which is consistent with the total dissociation energy of molecular 



2mH 



5.7 X 10^"^ erg where Xdis — 7.17 x 10 erg is the dissociation energy of H2 (ILiu et al.ll200l . 



This radius at the end of the simulation is about 2.5 times larger than the radius of the protostellar core obtained in 
[Masunaga & Inutsuka (200 0,) 4 j?0). and still increasi ng. The expansion in the adia batic accretion phase had be en 
already reported in iLarsonl (jl969D and also discussed in IStahler et al.l (jl986D (see also lWinkler fc NewmanI (|1980bD ). 
and our radius is consistent with their results. It continues to expand until almost all the gas in the first core 
has accreted onto the protostellar core when the accretion rate gets significantly low and the optical depth of the 
envelope becomes low enough for radiation cooling. It will take about the average free-fall time of the whole first core, 

^ff.FC = \J s.2g'p^c ~ 5.5 yrs where pre = = 1-5 x 10-i°gcm-3, Mpc ~ 3 x IO'^Mq and i?FC ~ 3 AU. Therefore 

this expansion is a highly transient phenomenon. In other words, our protostellar core has not settled yet. Note that 
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Figure 4. The evolution of the central gas density (top) and temperature (bottom) as functions of time. 



the gas does not flow outward in this expansion phase, but the newly accreted gas is loaded on top of the proto stellar 
core. It is different from the violent explosion, or "hiccup" of protostellar cores discussed years ago (e.g.. lBosslll989.) . 

5.2. Rotating Models 
5.2.1. Overview 

We show the evolution tracks of the central gas density and temperature as functions of time in Figure 01 and 
the evolution tracks in the p — T plane in Figure [5] The first cores are formed when the central density exceeds 

Pc ~ lO^^^gcm"^ and the second collapse begins when pc ~ 10~^gcm~^. The lifetimes of the first cores are about 
650, 720, 800, 850 and 950 years in SP, IS, IF, RS &nA RF, respectively. T he presence of rotation extends the first core 
lifetim e but its effect is not significant compared to non-magnetized cases (jSaigo et al.|[2Q08l: [Bati[20llt iTomida et al.l 
l2010al) . The resistive models have slightly longer lifetimes because magnetic fields are weakened by Ohmic dissipation 
and the efficiency of the angular momentum transport is reduced, but they are still not very long-lived. 

Interestingly, all the models show essentially the same thermal evolution (Figure [5]). This is because the central 
regions of the first cores contract almost spherically because of the efficient angul ar momenturn transport via m agnetic 
fields and their evolution converges on the Larson- Pension self-similar solution ()Larsonl [19691 : lPenstonlll9690 . In the 
central regions of the collapsing cloud cores, the thermal energy is dominantly larger than the magnetic energy and 
dissipation of magnetic fields (Joule heating) does not affect the thermal evolution significantly. Hence the thermal 
properties of the formed protostellar cores such as the gas entropy have already been set in the first core phase to 
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Figure 5. The evolution traeks of the central gas elements in the p — T plane. 



some extent and seem to depend weakly on some initial cloud parameters such as rotation and magnetic fields when 
there present efficient mechanisms of angular momentum transport. Without magnetic fields, the initial cond itions 
have larger impact on the thermal evolutions, particularly when fragmentation occurs due to fast initial rotation (jBatd 
1201 Ih . The thermal evolution may also depend on many oth er factors such as initial temperature, dust opacities and 
mass of the natal molecular cloud core (|Tomida et al.ll2010aD . 

5.2.2. Outflows and First Cores 

After the formation of the first c ores, slow and looscly-collimated outflows are launched from the first cores by 
magneto-centrifugal force (Blandford fc Payn e 1982; Kudoh ct al. 1998). We show the density and temperature dis- 
tributions of the outflow scale (corresponding to the grid level ^ = 8 or ^ 140 AU) and of the first core scale {I — 11 
or ~ 18 AU) in Figures [6] - [13] at the end of the first core phase. The outflows can be seen as slowly outgoing gas 
denser than the envelope that extends to z = 55, 70, 60 and 80 AU in IS, IF, RS and RF, respectively. The outflows 
are not prominent in the temperature plots except for the thin shock-heated layers because they are optically thin and 
radiation cooling/heating is efficient. We also show the profiles of the gas density, temperature and velocities along x- 
and z-axes in Figur e [H] Roughl v speaki ng, the properties of t he first cores and outflows are consistent with those in 
previous studies fCommercon et al. 2010: iTbmida et al]|2010bl etc.). 

The properties of the outflows such as velocities and traveling distances are similar in all the models. This is because 
the driving radii are similar in all the models, ~ 10 AU, where the gas density is not high enough for the resistivity 
to work. The outflow velocities are comparable to the rotational velocities at this radius, ^ Ikms^^. Therefore the 
traveling distances are almost proportional to the lifetimes of the first cores. It seems difficult to find the effect of the 
resistivity on the outflow pr operties, because the re is no essential difference in the outflows between the resistive and 
ideal MHD models (see also lYamada et al.ll2009l ). 

Although the outflows are similar, the first cores look quite different between the resistive and ideal MHD models. 
In the ideal models, the first cores are virtually non-rotating because of the efficient magnetic braking (i? ^ 3 AU in 
Figure [H)) . However, the magnetic braking is less efficient in the resistive models and there remains considerable 
amount of angular momentum. To show the effects of the resistivity clearly, we plot the total angular momentum 
within the first cores in Figure 1151 Here we simply measure the angular momentum where the gas density is above a 
critical value, pcrit = 10"^'^ g cm^"^. The resistivity is efficient where 10~^°gcm~'^. The magnetic flux is extracted 
from the central high density region and redistributed to the outer thinner region, but the redistributed magnetic fields 
do not show any significant effect dynamically at least in this phase. Both in the models with fast and slow rotation, 
the first cores in the resistive models attain about twice larger angular momenta compared with the corresponding 
ideal MHD models. Even these differences do not have significant impact on the evolution of the first cores, they 
become important later in the protostellar core phase. 

The outstanding difference between the ideal and resistive models is that the first cores and the surrounding (pseudo) 
disks warp in the ideal MHD models. These are likely to be the magnetic interchange instability (jSpruit fc Taam|[l990l : 
ISijruit et al..,1995; .Stehle fc Spruit .2001. : .Krasnopolskv et al..,2012. ). Uniformly-rotating disks are unstable when the 
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Figure 6. The vertical (top) and horizontal (bottom) cross sections of the gas density (left) and temperature (right) in the outflow scale 
(Z = 8 or ~ 140 AU) of Model 7F just before the second collapse. Projected velocity vectors are overplotted. 
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Figure 7. The same as Figure [6l but of IS. 
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Figure 8. The same as Figure [6l but of RF. 
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Figure 9. The same as Figure [8l but of RS. 

ratio between the poloidal magnetic flux density and the surface density decreases radially; '^'"^^J^'' < 0. Both ideal 

and resistive models satisfy this condition in 1 AU ^ i? 10 AU (Figure [TO)) , though only the ideal MHD models suffer 
from the warping in our simulations. There are several possibilities to explain this: the timescale of growth is long 
in the resistive cases or we need to consider more realistic situations including the effects of finite disk thickness and 
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Figure 11. The same as Figure [TOl but of IS. 
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Figure 14. The radial profiles of the gas density, temperature along the x- (in the disk mid-plane, red) and 2-axes (along the rotational 
axis, green), and the infall (red) and rotation (green) velocities along the x-axis (from top to bottom) at the end of the first core phase. 
From left to right, the different columns arc for Model IF, IS, RF and RS. 




Figure 15. The evolution of the total angular momenta in the first cores versus the central gas density. Magnetic fields are decoupled 
from fluid where Rm < 1 (yellow shaded region). 
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Figure 16. The "fiux-to-mass" ratio B^/T, (solid) and the infall velocity —Vr (dotted) of /^(red) and RF (blue) along the axis when 
the central density is pc ~ 4.1gcm~'', before the warping starts to grow. S and Bz is calculated in — 5AU < z < 5AU and mass- 
weighted average is used to evaluate Bz in this region. Both models satisfy the condition for the interchange instability, ''^^j/^^ < 0, in 

lAU£i?£ lOAU. This fi gure clearly shows that the magnetic fields are transported outward via Ohmic dissipation in the resistive MHD 
model. 

gas infall to understand these phenomena, but they are beyond the scope of this paper. Possibly, the magnetically 
driven warping instability (Lai 2003) may also contribute to the warping. We should note that even small vertical 
perturbation induces artificial reconnection of magnetic fields at the mid-plane because the magnetic fields are pinched 
onto the disk mid-plane and the directions of the magnetic fields above and below the mid-plane are anti-parallel, then 
things go chaotic. Because this warping does not affect the evolution of the protostellar cores, we do not discuss it 
further in this work. 

We can see that the gas within the dust evaporation front slowly infalls due to loss of the pressure gradient in all 
the models. However, the dust evaporation front is still confined in the first core and therefore its effects do not seem 
significa nt. This is different from the n on- magnetized RHD simulations where the front expands beyond the first core 
surface (jSchonke Sz Tscharnuterll2011| ). In our simulations, the angular momentum transport is very efficient and the 
first core properties are quite similar to the spherically symmetric cases, even in the resistive models. 

Interestingly, the size (height) of the first core is slightly larger in the resistive models. It is about 3 AU in the ideal 
MHD models, which is similar to the spherical model, but is about 5 AU in the resistive models. This is interpreted as 
a consequence of energy transport and additional heating by Ohmic dissipation. The magnetic fields are transported 
outward, then heat up and infiate the outer region (note that the resistivity is most effective around p ^ 10~^ gcm~'^). 
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Figure 17. The vertical (top) and horizontal (bottom) cross sections of the gas density (left) and temperature (right) in the protostellar 
core scale (I = 16 or ~ 0.54 AU) of Model IF at the end of the simulation. 
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Figure 18. The same as Figure [TTl but of IS. 
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Figure 20. The same as Figure [19] but of RS. 
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Figure 21. The same as Figure [141 but at the end of the simulations. 



5.2.3. Protostellar Cores and Jets 

We show the density and temperature cross sections of the the protostellar core scale (/ — 15 or 16, corresponding 
to ~ 1.1 AU or ^ 0.54 AU) at the end of the simulations in Figures [T71 - ?IU[ We also show the profiles along x- 
and z-axes in Figure [211 We stop our simulations when the central temperature reaches Tc ~ 10^ K, corresponding 
to 1.05,0.44,0.90 and 1.25 years after the formation of the protostellar cores in IF, IS, i?F and RS, respectively. Our 
simulations correspond to the earliest phases of the protostars. 

In the ideal MHD models, the properties of the protostellar cores are very similar to the spherical model. The 
evolution of the masses and radii of the protostellar cores are essentially identical to those in the spherical model SP. 
This is because the magnetic fields take almost all the angular momentum away from the gas in the central region. 
Contrarily, the resistive models have significantly larger angular momenta and the protostellar cores are strongly 
supported by rotation. To clarify the differences of the evolution between the models, we show the evolution of the 
radii, masses and angular momenta of the protostellar cores in Figure [2H Here we define the radius of the protostellar 
core as the radius where the infall velocity is largest (corresponding to the shock at the surface of the protostellar 
core) in the xy-plane, and measure the mass and angular momentum within this radius. The angular momenta in the 
protostellar cores in the resistive models are larger than those in the ideal MHD models by more than two orders of 
magnitudes. The rotationally-supported protostellar cores (or "circumstellar disks") are quickly built up within ^ 1 
year after the formation of the protostellar cores. At the end of the simulations, the radii of the disks are about 0.35 
AU in RF and 0.2 AU in RS, and they are still growing as the gas with larger angular momentum accretes. 

The protostellar cores in the resistive models also look like nearly spherical, but this is actually just a coincidence. 
They are supported by rotation in the horizontal direction, but they are vertically inflated due to the outflows. The 
toroidal magnetic fields are rapidly amplified in these rotating cores and the magnetic pressure gradient force drives the 
well-collimated outflows, or "jets" (Figures 1231 -[25]). The outflows are visible in the cross sections of the temperature 
as the hot towers. In the fast rotating model RF, the maximum outflow velocity reaches Vz ''^ 15kms~^, while it 
is ~ 6kms~^ in RS. These velocities are comparable to the rotational velocities seen in the protostellar cores, and 
therefore far faster than the outflows driven from the first cores because of the deeper gravitational potential. 

The magnetic fields in the rotating protostellar cores are quickly wound up and form the so-called magnetic tower. 
The tightly-wound magnetic fields are susceptible to the kink instability in long wavelengths. In our resistive models, 
the kink instability grows rapidly and the outflows start precession (Figures [23] - [^5t . Although this instability 
disturbs the coherent toroidal magnetic fields, the outfiow velocity is still getting accelerated because the bulk angular 
momentum in the protostellar core is increasing (Figure [22]) as the matter with higher angular momentum continuously 
accretes from the envelope, or the remnant of the first core. Therefore we expect that the outflow will be faster as the 
disk acquires larger angular momentum and the gravitational potential becomes deeper. 

In all the models, radiation feedback from the protostellar core formation on the first core and outer envelope is 
not significant. The first cores seem to remain almost unaffected after the formation of the protostellar cores. This is 
different from previous RHD simulations without magnetic fields; the surrounding first cores are heated up by radiation 
and bipolar outfiows are launched (|Bat&, 2010. .2011: .Schonke &: Tscharnutcr.2011, ). We need to study longer evolution 




Figure 22. The evolution of tlie radii, masses and angular momenta of the protostellar cores as functions of the central gas density. 
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Figure 23. 3D view of tlie protostellar core (I = 17) in Model RF, just after the formation of tlie protostellar core before the growth 
of the kink instability. The edge of the figure is ~ 0.27 AU. The left and bottom panels are cross sections of the density and the right 
panel shows the temperature cross section. The high density region (p > 10~^ gcm~^) is visualized with the orange surface. White arrows 
denote the direction of the fluid motion and white lines the magnetic field lines. Fast outflowing gas (v^ > 3kms~^) is volume-rendered 
with pale yellow. 




Figure 24. 3D view of the protostellar core (I = 16) in Model RF, in the growing phase of the kink instability. The edge of the figure is 
~ 0.54 AU. The gas with Vz > 4kms~^ is rendered with pale yellow. 
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Figure 25. 3D view of tlie protostellar core (I = 15) in Model RF at the end of the simulation, about three months after Figure [23l The 
kink instability is already grown significantly. The edge of the figure is ^ 1.1 AU. The gas with Vz > Tkms"'^ is rendered with pale yellow. 

to ad dress this problem because we only c alculate the earliest (i ^ 1 yr) evolution of the protostellar cores; iBatidMol 
1201 If ) and Sc honke fc Tscharnuterl ()2011[ ) followed the evolution of the protostellar cores more than 15 years, longer 
than the free-fall time of the first cores. 

5.3. Summary of Results 

We summarize the properties of the first cores and p r otoste llar cores at the end of the simulations in Table [2J The 
table also includes the results of iMasunaga fc Inutsukal ()200ClD and lMachida et all ()2008aD . The properties of the first 
cores and associated outflows in the rotating models are similar to those in the previous studies. The protostellar 
cores are significantly larger than those in the previous studies, but they are consequences of the transient expansion 
which happens in the earliest evolution of the protostellar cores. This transient expansion cannot be captured with 
the barotropic approximation in which the shock heating is not taken into account. The protostellar cores acquire 
their masses very quickly in this phase 0.02 Mq in a year). B ecause this phenomenon has very short time scale 
(estimated to be ~ 5 years), our models are not inconsistent with IMasunaga fc Inutsukal ()2000D . and also consistent 
with the present-day case of lOmukai eFall (|20Tot ). This expansion will affect the properties of the protostellar core 
when the core settles down after this expansion. It may be critically important in the further evolution of the protostars 
as the initial condition. However, to confidently discuss this phenomenon and the properties of the protostellar cores, 
we have to calculate the evolution far longer. 

The protostellar cores formed in the ideal MHD models are essentially non-rotating because of the efficient angular 
momentum transport via magnetic fields. In the resistive models, on the other hand, the protostellar cores attain 
considerably large angular momenta and the rotationally-supported disks emerge in their earliest phases, and they will 
continuously evolve into circumstellar disks. Although it is difficult to distinguish the disk from the central protostar 
in our simulations at this phase, they are going to separate into a thermally-supported protostar and a rotationally- 
supported disk when radiation cooling takes place and reduces the thermal support. Thus Ohmic dissipation remedies 
the magnetic braking catastrophe. 

Despite the significant difference of the thermal evolution in the protostellar cores, the properties of the outflows 
and circumstellar disks associated with the protos tellar cores in the resi stive models are consistent with the previous 
MHD studies using the barotropic approximation ('M achida et al.|[2008aD . This is because they are mainly determined 
by the interaction between magnetic fields and rotation, and the effects of the thermal evolution have less significant 
impacts on their properties. 

6. CONCLUSIONS AND DISCUSSIONS 

We performed 3D nested-grid RMHD simulations of formation of protostellar cores from molecular cloud cores with 
and without Ohmic dissipation, and revealed the earliest evolution (only 1 year after the formation) of the protostellar 
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Model 


Tpc (yrs) 


iJFc(AU) 


Dof(AU) 


-Rpc(Rg) 


Mpc(M0) 


Vjet (kms-1) 


SP 


650 


3 





10* 


0.02 





IS 


720 


3 


55 


10* 


0.02 





IF 


800 


3 


70 


17* 


0.02 





RS 


850 


5 


60 


45 


0.02 


5 


RF 


950 


5 


80 


75 


0.02 


15 


MI 


650 


3 





4 


0.016 





MR'i 




0.5 




8.2 


0.008 


15 



Table 2 

Summary of the properties of the first cores and the protostellar cores. Prom left to right: the lifetime of the first core, the 
thermally-supported radius of the first core, the distance the outflow traveled during the first core lifetime, the thermally- or 
rotationally-supportcd radius of the protostellar core, the mass of the protostellar core and the maximum velocity of the jet driven from 
the protostellar core. The quantities marked with * are strongly affected by the transient expansion. For comparison, we also show the 
results of Masunaga & Inutsuka (2000) {MI), the snapshot labeled "10" which corresponds to 19 years after the second collapse, and the 
model MR of Machida Ot al. ( 2008a ). f: We adopt their disk radius as the protostellar core radius and we show the total mass of the 
protostar, disk and outflows because we do not distinguish them in this study. The size of the thermally supported protostellar core in 

MR is small, ~ 1^0, because of the barotropic approximation. 

cores. These simulations are, to our knowledge, the first 3D RMHD simulations in the world following the whole 
evolution from molecular cloud cores to protostellar cores with realistic physics. We successfully revealed the realistic 
evolution in the earl y phase of protostellar c ollapse. 

Preceding studies (|Tomida et al.l l2010bl lat iBatd 120111) had showed that the barotropic approximation fails to re- 
produce the realistic thermal properties, and here we also revealed that the discrepancy is more prominent in the 
protostellar core phase, because the barotropic approximation does not take account of the shock heating and there- 
fore tend to underestimate the temperature, which results in the smaller radii of the first and protostellar cores and 
cannot reproduce the transient expansion of the protostellar cores. In order to capture the realistic thermal evolution, 
RMHD simulations are thus essential. Our result s can be used as th e initial and boundary conditions for the stellar 
evolution simulations (e.g., [Baraffc ct al. 2009; H osokawa et al.ir2011|) . 

We found two qualitatively different phenomena between the ideal MHD models and the resistive MHD models: 
formation of circumstellar disks and outflows from protostellar cores. In the resistive models, the circumstellar disks are 
quickly built up in the vicinity of the formed protostars. Although the formed disks are still ver y small (R < 0.35 AU) , 
we expect the disks will grow continuously as the gas with larger angular momentum accretes. iMachida et al.l (|2011l ) 
performed resistive MHD simulations with the sink cell technique and the barotropic approximation and showed 
that the rotationally-supported disk grows to R > 200 AU while the high density region where Ohmic dissipation 
works expands as the disk grows. The two-component outflows are launched from different scales as the natural by- 
products: the slow loosely-coUimated outflow driven from the flrst core and the fast well-coUimat ed jet driven from the 
protostellar core. These are consistent with previous studies using the barotropic approximation (jMachida et al.l[2006al . 
[2007l[2008ah . Although the maximum velocity of the outflow is still not very fast {v^ ~ 15kms ^), we expect that it 
will get faster as the protostellar core grows. Meanwhile, we expect that the outflow from the flrst core is continuously 
driven from the out er disk, or the remnant of the first core until the accretion stops (Machida & Matsumoto 201 IJ; 
IMachida et"al]l2011 1. If it is the case, our si mulations can na t urally explain the obse rved high-velocity jets, especially 
the multi-component outflows (for example. iLee et ani2000t iSantiago-Garcia et al.l[2009 'l. In the ideal MHD cases, 
on the other hand, the protostellar cores are almost the same with the spherically symmetric case and no outflow is 
launched as a result of efficient magnetic braking. These results do not necessarily mean that circumstellar disks are 
never formed in ideal MHD simulations. The circumstellar disks and outflows may be formed later when the gas with 
large angular momentum accretes sufficiently, because the magnetic braking is not a process which reduces the total 
angular momentum but transports the angular momentum from the disk to the thin envelope within the molecular 
cloud via Alfven waves. To simulate such a system, long-term simulations with reasonable boundary conditions are 
required. 

Although "when and how circumstellar disks are formed?" is a crucial question, there has been no conclusive answer. 
Our results suggest that a circumstellar disk is formed in the earliest phase of protostar evolution, essentially in parallel, 
[jjargcnscn ct al. (2009) suggested that Class-0 sources (young, embedded) are associated with more massive disks than 
Class-I sources, which support the early formation of circumstellar disks. We need more elaborate observations with 
higher spatial resolution to answer this question, and we expect the Atacama Large Millimeter/submillimeter Array 
(ALMA) will provide good clues to this problem. 

In our simulations, the formation of the protostellar core does not affect the first core and outer envelope. Considering 
the large energy r elease d in t he seco nd collapse and subsequ e nt acc retion with the high accretion rate, we expect that 
the feedback like IBatd (|201 1) and Schonke &: Tscharnuterl (|2011l ) reported may happen; bipolar outflows may be 
launched by radiation heating from the protostellar core. There are two possible interpretations of the discrepancy: 
magnetic fields and duration of the simulations. Because of the efficient angular momentum transport due to magnetic 
fields, the first cores are not supported by rotation even in the resistive models and the st ructu r es of the accretion 
flows are completely different between the models with and without magnetic fields. Indeed, IBatd ()2011f ) showed that 
the outflows are less energetic in the models with slower rotation. Therefore magnetic flelds possibly can prevent 
formation of such outflows driven by radiation heating from the protostellar cores. 
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The other point is that we o nly calculate about 1 year after the protostellar core formation while iBatd ()201lD and 
iSchonke fc Tscharnutej (|2011[ ) followed the evolution longer than 15 years. Because 1 year is shorter than the free-fall 
time of the first core, the protostellar cores are still deeply obscured in the remnants of the first cores in our simulations. 
If this is the case, how long will it take for the feedback to become prominent? As we mentioned above, the answer is 
the free-fall timescale of the first core, which is about 5 years. Another interesting factor is the interaction between 
the first core and the outflow from the protostellar core. When the outflow from the protostellar core penetrates 
the first core, it will have some influence on the first core. Assuming the constant traveling velocity of 15kms^^, it 
will reach the surface of the first core in about 1.5 years. Unfortunately, it is almost hopeless to follow such a long 
evolution with our current code because the simulation timesteps are too small, less than 1 minute, and are still getting 
smaller. In order to pe rform th e long-term simulations, we have to modify our code, for example, introducing the sink 
particle technique (Bat e et"an [T995: Krum holz ct al. 2 004; Fcdcrrath ct al. 201 0) to replace the protostellar core with 
a subgrid stellar evolution model (e.g., Baraffe et al]l2002i : lYorke fc Bodenheinieill2008t iHosokawa fc Omukaill2009f ). 
To construct reasonable subgrid models, our results can be used as templates. Since the radi ation feedback from the 
formed protostar c an be highly anisotropic ("flashlight effect". lYorke fc Bodenheimeii [19991 ) due to the small-scale 
optically-thick disk (jVaidva et al.|[2009t iTanaka fc Nakamotcill201lD . high resolution down to the protostellar core scale 
is of crucial importance even when we introduce the sink particle and subgrid stellar evolution model. 

All the physical processes involved in this work play crucial roles in star formation processes. Magnetic fields 
efhciently extract the angular momentum and drive the outflows from the first cores. Ohmic dissipation suppresses 
the angular momentum transport which is too efficient in the ideal MHD models and enables the formation of the 
circumstellar disks and the fast outflows from the protostellar cores. The EOS and radiation transfer give the realistic 
thermal evolution and enable us to discuss the dynamics and thermodynamics quantitatively, which cannot be achieved 
with the barotropic approximation. We should note, however, that there is no established test for such a complicated 
system, although we tested each part of our code separately with standard tests. To verify the validity of the codes, 
we are planning to perform comparisons with other groups. 
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APPENDIX 1. EQUATION OF STATE 

In some previous studies, the simple EOS of perfect gas has been often used, i.e., the adiabatic index F is assumed 
to be a constant throughout simulations, 7/5 which is valid for completely idealized diatomic molecular gas, or 5/3 
which is correct only when the gas temperature is very low. However, the adiabatic index is not constant in reality. In 
the low temperature molecular cloud (T ^ 100 K), molecular hydrogen behaves like monoatomic gas, F ~ 5/3, because 
the coUisional energy is insufficient to excite the rotational degrees of freedom and only translations are excited. When 
the gas temperature exceeds T ^ 100 K, the rotational degrees start to be excited and contribute to the heat capacity. 
Then the adiabatic index decreases, F ~ 7/5, close to that of the diatomic gas. The adiabatic index is of critical 
importance in the thermal evolut ion of the gas, and also in the stability of the gas against the gravitational instability 
(jStamatellos fc Whitworthll2009| ): the stiffer gas (with larger F) is more stable gravitationally because it reacts more 
strongly against compression. Therefore we need to consider these quantum thermodynamic behaviors in our EOS. 

The second collapse is driven by the endothermic reaction of hydrogen molecule dissociation. In order to simulate the 
evolution in the second collapse phase, we need to take the chemical reactions into account. However, it requires quite 
a large computational cost to solve the network of chemical reactions in every fluid element while solving (radiation) 
hydrodynamics. Fortunately, because we are mainly interested in dense gas, we can assume that the timescale of 
chemical reactions is shorter than the dynamical timescale. Therefore we implement the chemical reactions related to 
major species within the EOS on the assumption of the local thermodynamic and chemical equilibrium. 

We require some assumptions to calculate the EOS for simplicity: 

• The gas is in the local thermodynamic and chemical equilibrium (except for the ortho/para ratio of molecular 
hydrogen; see below). 

• All the atoms, molecules and ions are in the ground states. 

• The EOS can be calculated by simple summation of each component and each degree of freedom, i.e., the 
interactions between components and other non-ideal effects are neglected. 

• Only seven major species (H2, H, H+, He, He+, Hc^+ and e~) are considered and other elements are neglected. 

Based on these assumptions, we calculate the EOS using the statistical mechanics theory. Here we also assume that 
the gas has the solar abundance, X = 0.7 and Y — 0.28. 
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PARTITION FUNCTIONS 

Here we describe partition functions of each element. The partition function of a species i can be divided into 
five parts; translation Ztr,i, rotation Z^otA, vibration Z^ihA, spin Zgpin.i and electron excitation Zcicc,i (we include 
contributions from H2 dissociation and ionization of hydrogen and helium in the electron excitation part) : 



V X Ztr.t X Zi-c 



-'spin.z ^c\cc,i-! 



(45) 



where V is the volume. We calculate these partition functions by the standard procedure. In the following descriptions, 
rrii denotes the mass of i-species, h the Planck constant and k the Boltzmann constant, respectively. The partition 
functions not explicitly described are unity. 



Molecular hydrogen: 
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where ^rot = 170.64 K is the excitation temperature of rotation and 9vih = 5984.48 K is that of vibration. Molecular 
hydrogen has two forms: orthohydrogen with aligned nuclear spins and odd rotational states, and parahydrogen with 
antiparallel nuclear spins and even rotational states. Here we assumed the ratio of orthohydrogen to parahydrogen 
is 3:1. This ratio has significant impact on the dynamics of collapsing molecular cloud cores in the relatively low 
temperature region because thermodynamic properties related to rotation of molecular hydrogen depend on the nuclear 
spins and rotational states. Unfortunately this ratio in star forming regions is quite unclear because of observational 
difficulties. Although parahydrogen is more stable, some observations of int erstellar dark clouds suggest that the ratio 
is considerably far from the equilibrium value even i n the cold environme nt; iPagani et all ()2011|) proposed that the or- 
tho/para ratio is larger than 0.1. On the other hand. iDislaire et"al] ()2011[ ) claimed that the ratio is quite small, ~ 10"'^. 
In this work, we calculate the EOS using the ortho/para ratio of 3:1. This assum ption helps us interpret our simulation 
results because the adiabatic i ndex T decreases monotonically (Bolev et al. 2007) and also compare our results with re- 
cent simulations performed by lBatd p01Cli201L) (but iStamatellos fc Whitworth (,2009i ) assumed the equilibrium ratio). 



Atomic hydrogen: 



(2TnniikT) 



3/2 



^spin,H — 2 • - + 1 — 2, 



Ze,ec,H = 2exp(-|g), 

where Xdis = 7.17 x 10~^^ erg is the dissociation energy of H2 (|Liu et al.ll2009D . 
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Ionized hydrogen: 



^tr,H+ 



Z. 



spin,H'^ 



{2T:'mii+kT) 
2-i + l = 2, 



3/2 



Xdis + 2Yi, 



(56) 
(57) 



RMHD Simulations of Protostellar Core Formation 29 

where Xion — 2.18 x 10^^^ erg is the ionization energy of atomic hydrogen. 

Atomic and ionized hehum: ,„ ,^s-i/o 

„ (2^mHcfcr)3/2 

Ar,He — , loyj 

_ (27rmHe+fcr)3/2 

Ar,He+ . (oUj 

■^elec,He+ =exp i^^^^f-) ' (61) 

(2^mnc^iA:r)-V^ 

Ar,He2+ = p : (0^) 

^ / XHe,l + XHe,2\ .„„^ 
■^elec,He2+ =exp — . (63) 



where XHe.i = 3.94 x 10~^^ erg and XHe,2 = 8.72 x 10~^^ erg are the first and second ionization energies of hehum. 
Electron: 

Z.,,J^^^, (64) 



^spin,e = 2-- + l = 2. (65) 



CHEMICAL REACTIONS AND NUMBER DENSITIES 
The grand canonical partition function of each species is defined as: 

@,{HuV,T) = J^exp (^) ^ = exp kexp (m , (66) 



where /i, is the chemical potential of i-species and Ni is the total number of i-species. The grand potential can be 
derived from the grand canonical partition function: 

a(Mi, V,T) = -kTlog&i = -fcTZiexp (||) . (67) 

The total grand potential can be calculated from the summation of each component: 

fl = J2^i- (68) 

i 

We calculate required thermodynamic variables from these functions. First, we calculate the number density of each 
species based on the chemical equilibrium. The number density of i-species is derived from the partial derivative of 
the grand potential with respect to /ij: 

where Zi = Zi/V. This relation yields: 



7) ■ 

ft = fcTlog^. (70) 

Zi 

We consider only four reactions between the seven species dominant in relatively dense (but not too dense) gas for 
star formation problems: 

H2^2H, (71) 

H^H++c-, (72) 

He< — ^Hc++e", (73) 

Ee+< — ^Hc^+H-e". (74) 

Then the number densities can be calculated from the balance between the chemical potentials in these chemical 
reactions. 
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The RHS term of each equation, K^^ can be calculated from the partition functions. We have three additional relations; 
conservation of the total number density of hydrogen, conservation of the total number density of helium and the charge 
neutrality: 



2nH2 



-He 



We eliminate ?^H2 j "'H+ , "■Hc+ and nHc2+ from ((79l - [8T|) using ((75l - [78)) : 



TOHc 



(79) 

(80) 
(81) 



^ -?^Hc,l 
Ky{c,1 



-f^He,l-?^He,2 



^^lon — "tot 



(82) 
(83) 
(84) 



By substituting njj and njjc (we can determine the solution of 
to (|M)) . we obtain one non-linear equation of Ue'- 



uniquely since all the physical variables are positive) 
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Then we solve this equation numerically using the bi-section method (note that the first term of /(^e) is already 
modified to avoid the round-off error). Using the obtained rig, it is straightforward to calculate the number density of 
each species. 

In order to derive thermodynamic variables, the temperature and density derivatives of the number densities are 
required. We take logarithmic differentiation of ([75] -[5T |) . then they yield: 
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From this matrix equation we can numerically derive the required derivatives such as 
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The temperature derivatives of the partition functions can be calculated analytically. 
The total number density is 



n 



and its derivatives are 



9 Inn 
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dlnT 
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d\np 



THERMODYNAMIC VARIABLES 



We use the relation valid in ideal gas: 



P = nkT 



-kT 



where ji = is the mean molecular weight. The derivatives of the pressure are: 



ainP 
ainT 

91nP' 
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Pt = 1 + nr, 



Pp = np. 



The specific entropy of each species is derived from the grand potential: 



Si = - 



1 



and its derivatives are: 



dT 

dp 



w 

krij 



rflnT2 

9 In Ui 



dhip 



1 



krii 
P 

9 In rii 
dhiT 

dlnzi 
rflnT 



dhiZi 
dlnT 



din ; 



JJj_ 
kT 



dlnT 



1 



n£ 
P 



dT pT' 



(89) 
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Because the last terms in these derivatives are canceled out by taking summation of species when the chemical 
are in equilibrium, we can omit these terms. Then the total entropy and its derivatives are defined as: 
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reactions 

(98) 
(99) 
(100) 



Now we can derive thermodynamic properties we require in radiation hydrodynamic simulations. 
The isothermal and adiabatic sound speeds: 
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— ) -\I-Pp--^- (102) 

The adiabatic index: 

d\np J g P 



The internal energy per volume and its derivative: 



Ed In Zi 



de \ 

-g^j ^Cv^pTSt. (105) 

Note that we do n ot use the relation e = CyT w hich is valid only in the completely ideal cases, as IBolev et all (12001 
suggested (see also iBlack fc Bodenheimed ()1975l )). 

We tabularize these thermodynamic variables as functions of {p,eg) and {p,T) with sufficiently high resolution 
(Alogp = 0.05, A log eg = 0.025 and AlogT = 0.02) in p = fO^^a _ lOgcm^^ and T = 3 - f 0^ K. We use this EOS 
table with bi-log-linear interpolation. 

COMMENTS ON EOS 

Our treatment of the EOS for hydrogen and helium is valid in relatively low-density regions ranging from interstellar 
gas to the second collapse phase. However, in very dense region like the deep interior of the protostellar core, non-ideal 
effects are not negligible: interactions between particles, weak quantum effects in low-temperature but high-density 
region, pressure ionization of hydrogen, and contributions from other chemical species. Such non-ideal effects will 
affect the thermodynamics and the dynamics (e.g., the quasi-equilibrium state of the second core may vary). Actually, 
our EOS results in a serious unphysical behavior in the very high density region (p > 0.1 gcm"'^) that almost all the 
hydrogen particles are turned into the molecular form even when the gas temperature is high enough to destruct the 
hydrogen molecule. This is because of the assumption of the ideal che mical equilibrium, bu t in reality this assumption 
broke down there due to the neglected interactions between particles (jSaumon "eFaH [19951) . Then our EOS gives the 
considerably soft adiabatic index F (Figure [T]) because of the contribution from the vibration transitions of molecular 
hydrogen. Si nce this behavior is co mpletely unphysical, our EOS is invalid in such a high density region. More realistic 
EOS such as iSaumon et al.l (|1995l ) is required to calculate this region properly. 

APPENDIX 2. OPACITIES 

For the gray radiatio n transfer, we use the compiled tables of the Rosseland and Planck mean opacities of 
iSemenov et all (pOOl FI. iFerguson etHI (pOOSf PI and the Opacity Project (OP) (jSeaton et al.l IT991 R. For dust 
opacities, we adopt the composite aggregate dust model of the normal abundance from S emenov's mean opacit y ta- 
bles. Though Seme nov's tables also contain the gas opacities, we adopt the gas opacities of IFerguson et all ()2005D and 
iSeaton et al.l (I1994D because Semenov's P lanck mean opacity seems to be significantly lower than other opacity tables 
([Ferguson et al.l 120051 : ISeaton et al.lll99"^ where molecular and atomic lines dominate the opacity sources. Therefore 
we connect the tables of Semenov et al. (2003) and Ferguson et al. (2005) at the temperature where all the dust 
components evaporate. The dust evaporation temperature (it weakly depends on the gas d ensity, but in the typ ical 
density region, T ~ 1400 - 1500 K) is given in Semenov's opacity calculation code. We use IFerguson et al.l (|2005fl in 
the low temperature region (logT < 4.5) and OP in the high temperature region (4.5 < logT < 6). We tabulate these 
tables as functions of {p, T) and use them with bi-log-linear interpolation. 

Unfortunately, the opacity tables do not cover the whole required region. Figure [26l shows the coverage of the opacity 
tables in the p — T plane. The typical evolution track of the central gas element and the profile in the spherically 
symmet ric collapse are also p lotted. The dust opacities of ISemenov et al.l ([2003t ) cover 10^^^ < p(gcm~^) < 10~^. 
OP and [Ferguson et al.l ([20051 1 cover -8 < logi? < 1 where R = p/T^ and Tg =T(K)/10^. It is not serious that the 
very low density region is not covered because that region is extremely optically thin and the details of dust opacities 
do not matter there. We simply extrapolate the opacities by taking the nearest value at the given temperature. The 
high density region is far more problematic; we do not have proper opacities for the protostellar core. But actually 
the thermal evolution in this region is dominated by chemical reactions (dissociation and ionization) and radiation 
transfer is of less importance there because the gas is extremely optically thick. Therefore we dare to extrapolate the 
tables in the same manneJ*^. 

We show the combined opacity tables in Figure[27|as functions of p and T. Note that there is quite large uncertainty in 
opacities, especially due to the dust models such as the structure, composition, size distribution and so on. The thermal 

® http: / / www.mpia.de/homes /henning/Dust_opacities / Opacities /opacities. html 
^ http:/ /webs. wichita.edu/physics/opacity/ 

http:/ /cdsweb. u-strasbg.fr/topbase/TheOP.htmI 

We must be careful, however, that this problem becomes more serious when we calculate the evolution of the protostar longer than 
the Kelvin- Helmholtz timescale before the ons et of convection (or d euterium burning). As we can see from Figure [251 the earliest phase 
of the protostar is convectively stable (see also lStahler et aTl 11980al fEDl and therefore radiation plays a critical role in heat transfer. More 
elaborate opacity tables covering wider region are highly demanded. 
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Figure 26. The blue shaded re gion is covered by ga s opacity tables llSeaton et al.in"994l : IFereuson et al.l 120051 ) and the yellow region 
is covered by dust opacity tables l lSemenov et al.ll2003ll . The border line between blue and yellow corresponds to the dust evaporation 
temperature. The red line represents the typical evolution track of the central gas element in the spherical protostellar collapse and the 
green line does the profile at the end of the simulation. 
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Figure 27. The Rosseland (left) and Planck (right) mean opacities. 

evolution and dynamics in star formation remain qualitatively similar even when we change the dust parameters, but 
the observational properties such as Spectral Energy Distributions (SED) are directly affected by the differences 
between the (monochromatic) opacities. 

APPENDIX 3. RESISTIVITY 

We calculate Ohmic resistivity considering both thermal and non-thermal particles. For non-thermal processes, we 
adopt the table of resistivity based on the reaction network model between dust grains and chemical species constructed 
by Umebayashi & Nakano (2009) (see also Okuzvmri (2009)). The table of tint is given as a function of the gas density 
p, the temperature T and the ionization rate ^. Here we assume the typical interstellar ionization rate due to the 
cosmic rays, £,CR ~ 10~"'^^s~^, neglecting shielding by the gas. Since the attenuation depth of the cosmic rays is about 
lOOgcm"^, the ionization rate will be lower in the deep interior of the first core. In this sense, and also because 
we neglect ambipolar diffusion, our models corresponds to the lower limit of (but still highly efficient) magnetic flux 
loss. Because decay of radionuclides (^^Al is the dominant source) considerably contributes to the ionization rate. 
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Figure 28. The resistivity rj and the magnetic Reynolds number Rr, 
decoupled from fluid where Rm < 1 (yellow shaded region). 



are plotted as functions of the gas density. Magnetic fields are 



£,RA ~ 7.6 X 10"^^ ()Uniebavashi fc Nakanci|[2009[ ). the effect from neglecting the shielding is as large as about an 
order of magnitude at most. We should mention that there is still large uncertainty in the resistivity from the grain 
properties such as the structure, composition, size distribution and so on. 

In order to calculate the resistivity for our simulations, it is sufficient to estimate the contribution to the thermal 
ionization from the species which has low ionization energy. Potassium (K) has very low ionization energy {kT^Q^ ^ 
4.33 eV) and sufhciently abundant to recover good coupling between gas and magnetic fields, therefore it is the most 
important electron-supplying species in star formation processes. Here we calculate the resistivity due to the thermal 
ionization of potassium on the assumption of thermal equilibrium using the following equation: 



?7t = 7.5 X 10 exp 



/25188K 



r-i/y/' cm^s-i. 



(106) 



We calculate the total resistivity as follows: 



Vnt' 



(107) 



because the resistivity is inversely proportional to the ionization degree. 
We show the resistivity and magnetic Reynolds number Rm = vgXj/r] as a function of the gas density in Figure 

■^'^ and the free-fall velocity vg 



Here we consider the Jeans length Aj = 27rc< 



- /GMj 



(where Mj 



32Gp "-^'^ ""^ y^.^^.^j — y - g-Ajp) 

as typical length and velocity scales to estimate the magnetic Reynolds number. To draw this plot, we adopt the 
barotropic approximation as a typical thermal evolution to calculate the gas temperature: 

r-H 

T = max 



10,10 X 



(- 

\ Pcrit 



K, 



(108) 



where pcrit = 2 x 10 ^'^gcm ^ is the critical density and F = 7/5 is the adiabatic index. The resistivity steeply 
decreases in p ^ 10~* g cm~'^ because of the thermal ionization of potassium. The magnetic fields are decoupled from 
fluid where the magnetic Re ynolds nu mber is less than unity. Our r e sistivity is significa ntly lower than the resistivity of 
iKunz &: MouschoviasI ()2009[ ) (see also lKunz fc MouschoviasI ()2010D : iDapp et all ()2012|) ). This difference mainly comes 
from the ionization rates; they assume that decay of '"'K (^k ~ 2.43 x 10"^'^ s~^) is the dominant ionization source in 
the dense region where cosmic rays cannot pe netrate, and neglect th e contribution from short-lived species like ^^Al. 
Our resistivity is quite similar to that used in lMachida et al.l ( 2006al ). 
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